Spatial Regression

Jacob Patterson-Stein

Agenda

Key points
Does space matter? A quick simulation
A little bit (more) on regression
A real example
Resources
Questions

Key Points of This Presentation

If you have spatial coordinates or you think there is a theoretical reason why results might be spatially related, you should conduct spatial correlation tests.
Data simulation is an important way to test model specification and understand what might be going on in your data.
Always visualize.

Estimating a relationship

Let’s imagine we have a treatment and some outcome, say coffee (treatment) and work performance (outcome).
Literature, prior experience, and our boss suggests this should be a strongly positive, statistically significant relationship.
We fit a model to the data to estimate the relationship.
We get a result!

So you get some data

When things don’t work how they should

Let’s say you do all that and look at your model output and the estimate is…not significant!

Why we might have non-significant results

Could be an issue with the number of observations
Could be that we have poor model fit
Could be omitted variable bias
Indeed, there could be any number of things that can plague any model

Simulation is an important tool for model testing

To quote Gelman, et al (2020), simulation helps you understand if you have the right model under different scenarios.
It also allows you to trouble shoot and better understand what is going on.
You can create a “true” effect and then see how well your model actually identifies this effect. If it can’t identify this effect, you have more work to do.

Set some criteria around our question

Specify our model \[ \operatorname{Work} = \alpha + \beta_{1}(\operatorname{Coffee}) + \epsilon \] Where, the average effect of coffee, \(beta_1\), is set at 0.1, and \(\epsilon\) is a random error, normally distributed with mean 0, sd of 1. If our model is correctly specified, we should recover this effect.

Run the model a few thousand times

What does this mean?

Our model appears to be generally ok. So what else is going on?
Recall that we also have latitude and longitude columns in our data.
It doesn’t seem like we have high amounts of correlation across variables, so there might be something going on in the residual error.

Start with the eye test and move to the I test!

Let’s look at the residuals of our model along with the spatial lag of residuals (i.e., the residuals of nearest neighbors).

Moving to Moran’s I Test

Moran’s I test “measures spatial autocorrelation based on both feature locations and feature values simultaneously.” Basically, it is a measure of how similar each unit’s residual is with some \(k\) set of neighbors.

In math: \[ I = \frac{N}{W} \cdot \frac{\sum_{i=1}^{N} \sum_{j=1}^{N} w_{ij} (x_i - \bar{x})(x_j - \bar{x})}{\sum_{i=1}^{N} (x_i - \bar{x})^2} \] Where, \(N\) is the sample size for unit \(i\) across area \(j\), \(x\) is the variable of interest, \(W\) is the weight created through row-wise standardization so that the sum of all of unit \(i\)’s neighbor’s weights is equal to 1. The ArcGIS documentation on this is very good.

Run Moran’s I

Let’s start with a test that assesses the correlation with a given unit’s 15 nearest neighbors. Moran’s I is based on a basic hypothesis testing framework, where:

\(H_0: I = E[I]\), a null of no spatial autocorrelation, i.e., I (the correlation between neighbors) is equal to the expected relationship, with, \(E[I] = \frac{-1}{n-1}\)
\(H_1: I \neq E[I]\), an alternative hypothesis that spatial correlation is present

A positive Moran’s I suggests positive clustering, i.e., units tend to have similar values, while a negative value suggests dispersion different from random.


    Global Moran I for regression residuals

data:  
model: lm(formula = Work ~ Coffee + Age + Income, data = data)
weights: listw

Moran I statistic standard deviate = 13.552, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
     0.422427205     -0.010046388      0.001018314

More on Moran’s I

When we talk about p-values, what we are talking about is the probability of obtaining a Moran’s I at least as extreme as the one observed assuming the null is true for our sample. The null implies random distribution of residuals, or the middle box below.

Let’s run our model but with a spatial error regression

Quick recap

The spatial error model is essentially adding a weighted residual to account for the fact that there are relationships not controlled for in your original specification.

\[ {\mathbf y} = {\mathbf X}{\mathbf \beta} + {\mathbf u}, \qquad {\mathbf u} = \rho_{\mathrm{Err}} {\mathbf W} {\mathbf u} + {\mathbf \varepsilon} \]

Put simply, sometimes the omitted variable biasing your results is right there next to you, and your neighbor, and your neighbor’s neighbor.

A real example

This is a map of depression in Seattle, Washington. We want to understand the relationship between neighborhood characteristics and depression prevalence.

Do we have spatial autocorrelation?

We can map the residuals of a basic linear model to get a better idea of correlation

Spatial autocorrelation in Seattle

Check Moran’s I


    Global Moran I for regression residuals

data:  
model: lm(formula = DEP_CrudePrev ~ unempr + pmob + pcol + ppov +
pnhblk + phisp + log(tpop), data = sea.tracts)
weights: seaw

Moran I statistic standard deviate = 6.6676, p-value = 1.3e-11
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
     0.309531240     -0.026084108      0.002533614

Regression results


<table style="text-align:center"><caption><strong>Title: Regression Results</strong></caption>
<tr><td colspan="3" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="2"><em>Dependent variable:</em></td></tr>
<tr><td></td><td colspan="2" style="border-bottom: 1px solid black"></td></tr>
<tr><td style="text-align:left"></td><td colspan="2">DEP_CrudePrev</td></tr>
<tr><td style="text-align:left"></td><td><em>OLS</em></td><td><em>spatial</em></td></tr>
<tr><td style="text-align:left"></td><td><em></em></td><td><em>error</em></td></tr>
<tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td></tr>
<tr><td colspan="3" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">unempr</td><td>-0.008</td><td>0.019</td></tr>
<tr><td style="text-align:left"></td><td>(0.071)</td><td>(0.048)</td></tr>
<tr><td style="text-align:left"></td><td></td><td></td></tr>
<tr><td style="text-align:left">pmob</td><td>0.006</td><td>0.047<sup>***</sup></td></tr>
<tr><td style="text-align:left"></td><td>(0.017)</td><td>(0.015)</td></tr>
<tr><td style="text-align:left"></td><td></td><td></td></tr>
<tr><td style="text-align:left">pcol</td><td>0.045<sup>***</sup></td><td>0.010</td></tr>
<tr><td style="text-align:left"></td><td>(0.014)</td><td>(0.014)</td></tr>
<tr><td style="text-align:left"></td><td></td><td></td></tr>
<tr><td style="text-align:left">ppov</td><td>0.117<sup>***</sup></td><td>0.081<sup>***</sup></td></tr>
<tr><td style="text-align:left"></td><td>(0.023)</td><td>(0.019)</td></tr>
<tr><td style="text-align:left"></td><td></td><td></td></tr>
<tr><td style="text-align:left">pnhblk</td><td>-0.083<sup>***</sup></td><td>0.006</td></tr>
<tr><td style="text-align:left"></td><td>(0.023)</td><td>(0.020)</td></tr>
<tr><td style="text-align:left"></td><td></td><td></td></tr>
<tr><td style="text-align:left">phisp</td><td>0.087<sup>**</sup></td><td>0.037</td></tr>
<tr><td style="text-align:left"></td><td>(0.034)</td><td>(0.026)</td></tr>
<tr><td style="text-align:left"></td><td></td><td></td></tr>
<tr><td style="text-align:left">log(tpop)</td><td>0.386</td><td>-0.440</td></tr>
<tr><td style="text-align:left"></td><td>(0.388)</td><td>(0.279)</td></tr>
<tr><td style="text-align:left"></td><td></td><td></td></tr>
<tr><td style="text-align:left">Constant</td><td>15.909<sup>***</sup></td><td>24.377<sup>***</sup></td></tr>
<tr><td style="text-align:left"></td><td>(3.247)</td><td>(2.519)</td></tr>
<tr><td style="text-align:left"></td><td></td><td></td></tr>
<tr><td colspan="3" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>134</td><td>134</td></tr>
<tr><td style="text-align:left">R<sup>2</sup></td><td>0.425</td><td></td></tr>
<tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.393</td><td></td></tr>
<tr><td style="text-align:left">Log Likelihood</td><td></td><td>-216.043</td></tr>
<tr><td style="text-align:left">sigma<sup>2</sup></td><td></td><td>1.233</td></tr>
<tr><td style="text-align:left">Akaike Inf. Crit.</td><td></td><td>452.087</td></tr>
<tr><td style="text-align:left">Residual Std. Error</td><td>1.560 (df = 126)</td><td></td></tr>
<tr><td style="text-align:left">F Statistic</td><td>13.284<sup>***</sup> (df = 7; 126)</td><td></td></tr>
<tr><td style="text-align:left">Wald Test</td><td></td><td>205.489<sup>***</sup> (df = 1)</td></tr>
<tr><td style="text-align:left">LR Test</td><td></td><td>59.094<sup>***</sup> (df = 1)</td></tr>
<tr><td colspan="3" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="2" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
</table>

More advanced stuff and other materials

Thank you!